Matlab
code 2.4: Matlab file “Figure 2-9.m”
%--------------------------------------------------------------------
% This code can be used to generate Figure 2.9
%--------------------------------------------------------------------
clear all;
close all;
%% the theoretical curve of rotating target
micro-Doppler characteristic when radar platform vibrating
c = 3e8;
j = sqrt(-1);
fc = 10e9; % carrier frequency of transmitted
signal
v = 0; % translational velocity of target
cord = 1000*[300 100 500]; % coordinates of
local coordinate system's origin in the radar coordinate system
colo = [0 0 0]; % coordinates of radar in the
radar coordinate system
R0 = cord-colo;
n = R0/sqrt(sum(R0.^2)); % unit vector of
Line-of-Sight(LOS)
ae = [0 pi/4 pi/5]; % Euler angle
tgt = [0 0 0;3 1.5 1.5;-3 -1.5 -1.5]; % three
point scatterers
w = [pi 2*pi pi]; % angular velocity of
rotation
T = 0.8165; % rotating period
ri1 = [cos(ae(3)) sin(ae(3)) 0;-sin(ae(3))
cos(ae(3)) 0;0 0 1];
ri2 = [1 0 0;0 cos(ae(2)) sin(ae(2));0
-sin(ae(2)) cos(ae(2))];
ri3 = [cos(ae(1)) sin(ae(1)) 0;-sin(ae(1))
cos(ae(1)) 0;0 0 1];
ri = ri1*ri2*ri3; % initial rotation matrix
tgtr = ri*tgt; % scatterers in the reference
coordinate system
alpha = atan(cord(2)/cord(1)); % azimuth angle
beita =
atan(cord(3)/(sqrt(cord(1)^2+cord(2)^2))); % elevation angle
alpha0 = pi/10; % azimuth angle of radar
platform vibrating axis
beita0 = pi/3; % elevation angle of radar
platform vibrating axis
dr = 0.02; % vibration amplitude of radar
platform
fr = 5; % vibration frequency of radar platform
tn = length(w); % number of scatterers
r0 = zeros(3);
for i = 1:tn
r0(i,:) = tgt(i,:)-colo;
end
r0r = ri*r0'; % scatterers in the reference
coordinate system
w = ri*w'; % angular velocity of rotation in
the reference coordinate system
omega = sqrt(sum(w.^2));
we = w/omega; % unit vector of rotation
wr = [0 -we(3) we(2);we(3) 0 -we(1);-we(2)
we(1) 0]; % skew symmetric matrix
t = 2; % radar illumimated time
prf = 3000; % pulse repetition frequency
pri = 1/prf; % pulse repetition interval
dt = 0:pri:t-pri; % time sampling interval
m = length(dt);
fmd1 = zeros(length(ae),m); % theoretical
micro-Doppler frequency
fmd2 = zeros(length(ae),m);
for i = 1:m
fmd1(:,i) =
(2*fc*omega/c)*((wr^2.*sin(omega*dt(i))+wr.*cos(omega*dt(i)))*ri*(r0'))'*n'...
-(4*pi*fc*fr*dr/c)*(cos(alpha-alpha0)*cos(beita)*cos(beita0)...
+sin(beita)*sin(beita0))*cos(2*pi*fr*dt(i));
fmd2(:,i) =
(2*fc*omega/c)*((wr^2.*sin(omega*dt(i))+wr.*cos(omega*dt(i)))*ri*(r0'))'*n';
end
figure(1)
plot(dt,fmd1,'r') % when radar platform
vibrating
hold on
plot(dt,fmd2,'b') % without radar platform
vibrating
xlabel('Time (s)')
ylabel('Frequency (Hz)')
axis([0,2,-1500,1500])
%% the time-frequency analysis result of
rotating target micro-Doppler characteristic when radar platform vibrating
r = zeros(length(tgt(:,1)),length(dt)); %
distance between the scatterers and radar
for i = 1:m
for n
= 1:length(tgt(:,1))
r(n,i) =
sqrt((R0'+v*t+((eye(3)+wr*sin(omega*dt(i))+wr^2*(1-cos(omega*dt(i))))*r0r(:,n))...
-(dr*sin(2*pi*fr*dt(i))*[cos(alpha0)*cos(beita0);sin(alpha0)*cos(beita0);sin(beita0)]))'...
*(R0'+v*t+((eye(3)+wr*sin(omega*dt(i))+wr^2*(1-cos(omega*dt(i))))*r0r(:,n))...
-(dr*sin(2*pi*fr*dt(i))*[cos(alpha0)*cos(beita0);sin(alpha0)*cos(beita0);sin(beita0)])));
end
end
s = zeros(length(dt),length(tgt(:,1))); % echo
signal matrix
for i = 1:length(tgt(:,1))
s(:,i) = exp(j*2*pi*fc*2*r(i,:)'/c);
end
st = sum(s,2);
N = 200; % number of Gabor coefficients in time
Q = 50; % degree of oversampling
tfr = tfrgabor(st.',N,Q); % Gabor transform
tfr_r = fftshift(tfr,1);
figure(2)
contour(linspace(min(dt),max(dt),length(tfr_r(1,:))),linspace(-prf/2,prf/2-1,length(tfr_r(:,1))),tfr_r,30)
xlabel('Time (s)')
ylabel('Frequency (Hz)')
axis([0,2,-1500,1500])